Bayesian GLMM Part5

1 Preparations

Load the necessary libraries

library(tidyverse) # for data wrangling
library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(knitr) # for kable
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(brms)
library(broom.mixed)
library(tidybayes)
library(bayesplot)
library(standist) # for visualizing distributions
library(rstanarm)
library(cmdstanr) # for cmdstan
library(ggeffects)
library(rstan)
library(DHARMa)
library(ggridges)
library(easystats) # framework for stats, modelling and visualisation
library(patchwork)
library(modelsummary)
source("helperFunctions.R")

2 Scenario

Some ornithologists were interested in the degree of sibling negotiations in owl chicks. Specifically, they wanted to explore how sibling negotiations were affected by feeding satiety and the sex of the parent returning to the nest. The ornithologists had accessed to a number of owl nests and were able to count (via recording equipment) the number of sibling negotiations (calls) that the owl chicks made when the parent returned to the nest.

We could hypothesise that the chicks might call more if they were hungry. As part of the investigation, the researchers were able to provided supplementary food. As such, they were able to manipulate the conditions such that sometimes the chicks in a nest would be considered deprived of supplementary food and at other times they were satiated.

As a parent returned, the researchers recorded the number of sibling negotiations (calls) along with the sex of the parent. Since the number of calls is likely to be a function of the number of chicks (the more chicks the more calls), the researchers also counted the number of siblings in the brood.

Each nest was measured on multiple occasions. Hence, we must include the nest as a random effect to account for the lack of independence between observations on the same set of siblings.

3 Read in the data

owls <- read_csv("../data/owls.csv", trim_ws = TRUE)

Rows: 599 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): Nest, FoodTreatment, SexParent
dbl (4): ArrivalTime, SiblingNegotiation, BroodSize, NegPerChick

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse

glimpse(owls)

Rows: 599
Columns: 7
$ Nest               <chr> "AutavauxTV", "AutavauxTV", "AutavauxTV", "Autavaux…
$ FoodTreatment      <chr> "Deprived", "Satiated", "Deprived", "Deprived", "De…
$ SexParent          <chr> "Male", "Male", "Male", "Male", "Male", "Male", "Ma…
$ ArrivalTime        <dbl> 22.25, 22.38, 22.53, 22.56, 22.61, 22.65, 22.76, 22…
$ SiblingNegotiation <dbl> 4, 0, 2, 2, 2, 2, 18, 4, 18, 0, 0, 3, 0, 3, 3, 6, 0…
$ BroodSize          <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
$ NegPerChick        <dbl> 0.8, 0.0, 0.4, 0.4, 0.4, 0.4, 3.6, 0.8, 3.6, 0.0, 0…

head

## Explore the first 6 rows of the data
head(owls)

# A tibble: 6 × 7
  Nest       FoodTreatment SexParent ArrivalTime SiblingNegotiation BroodSize
  <chr>      <chr>         <chr>           <dbl>              <dbl>     <dbl>
1 AutavauxTV Deprived      Male             22.2                  4         5
2 AutavauxTV Satiated      Male             22.4                  0         5
3 AutavauxTV Deprived      Male             22.5                  2         5
4 AutavauxTV Deprived      Male             22.6                  2         5
5 AutavauxTV Deprived      Male             22.6                  2         5
6 AutavauxTV Deprived      Male             22.6                  2         5
# ℹ 1 more variable: NegPerChick <dbl>

str

str(owls)

spc_tbl_ [599 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ Nest              : chr [1:599] "AutavauxTV" "AutavauxTV" "AutavauxTV" "AutavauxTV" ...
 $ FoodTreatment     : chr [1:599] "Deprived" "Satiated" "Deprived" "Deprived" ...
 $ SexParent         : chr [1:599] "Male" "Male" "Male" "Male" ...
 $ ArrivalTime       : num [1:599] 22.2 22.4 22.5 22.6 22.6 ...
 $ SiblingNegotiation: num [1:599] 4 0 2 2 2 2 18 4 18 0 ...
 $ BroodSize         : num [1:599] 5 5 5 5 5 5 5 5 5 5 ...
 $ NegPerChick       : num [1:599] 0.8 0 0.4 0.4 0.4 0.4 3.6 0.8 3.6 0 ...
 - attr(*, "spec")=
  .. cols(
  ..   Nest = col_character(),
  ..   FoodTreatment = col_character(),
  ..   SexParent = col_character(),
  ..   ArrivalTime = col_double(),
  ..   SiblingNegotiation = col_double(),
  ..   BroodSize = col_double(),
  ..   NegPerChick = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

Easystats (datawizard)

owls |> datawizard::data_codebook()

owls (599 rows and 7 variables, 7 shown)

ID | Name               | Type      | Missings |          Values |           N
---+--------------------+-----------+----------+-----------------+------------
1  | Nest               | character | 0 (0.0%) |      AutavauxTV |  28 ( 4.7%)
   |                    |           |          |          Bochet |  23 ( 3.8%)
   |                    |           |          |     Champmartin |  30 ( 5.0%)
   |                    |           |          |         ChEsard |  20 ( 3.3%)
   |                    |           |          |        Chevroux |  10 ( 1.7%)
   |                    |           |          | CorcellesFavres |  12 ( 2.0%)
   |                    |           |          |        Etrabloz |  34 ( 5.7%)
   |                    |           |          |           Forel |   4 ( 0.7%)
   |                    |           |          |          Franex |  26 ( 4.3%)
   |                    |           |          |            GDLV |  10 ( 1.7%)
   |                    |           |          |           (...) |            
---+--------------------+-----------+----------+-----------------+------------
2  | FoodTreatment      | character | 0 (0.0%) |        Deprived | 320 (53.4%)
   |                    |           |          |        Satiated | 279 (46.6%)
---+--------------------+-----------+----------+-----------------+------------
3  | SexParent          | character | 0 (0.0%) |          Female | 245 (40.9%)
   |                    |           |          |            Male | 354 (59.1%)
---+--------------------+-----------+----------+-----------------+------------
4  | ArrivalTime        | numeric   | 0 (0.0%) |  [21.71, 29.25] |         599
---+--------------------+-----------+----------+-----------------+------------
5  | SiblingNegotiation | numeric   | 0 (0.0%) |         [0, 32] |         599
---+--------------------+-----------+----------+-----------------+------------
6  | BroodSize          | numeric   | 0 (0.0%) |          [1, 7] |         599
---+--------------------+-----------+----------+-----------------+------------
7  | NegPerChick        | numeric   | 0 (0.0%) |        [0, 8.5] |         599
------------------------------------------------------------------------------

Skim (modelsummary)

owls |> modelsummary::datasummary_skim()

	Unique	Missing Pct.	Mean	SD	Min	Median	Max	Histogram
ArrivalTime	317	0	24.8	1.9	21.7	24.4	29.2
SiblingNegotiation	29	0	6.7	6.7	0.0	5.0	32.0
BroodSize	7	0	4.4	1.2	1.0	4.0	7.0
NegPerChick	78	0	1.6	1.6	0.0	1.2	8.5
		N	%
Nest	AutavauxTV	28	4.7
	Bochet	23	3.8
	Champmartin	30	5.0
	ChEsard	20	3.3
	Chevroux	10	1.7
	CorcellesFavres	12	2.0
	Etrabloz	34	5.7
	Forel	4	0.7
	Franex	26	4.3
	GDLV	10	1.7
	Gletterens	15	2.5
	Henniez	13	2.2
	Jeuss	19	3.2
	LesPlanches	17	2.8
	Lucens	29	4.8
	Lully	17	2.8
	Marnand	27	4.5
	Montet	41	6.8
	Murist	24	4.0
	Oleyes	52	8.7
	Payerne	25	4.2
	Rueyes	17	2.8
	Seiry	26	4.3
	Sevaz	4	0.7
	StAubin	23	3.8
	Trey	19	3.2
	Yvonnand	34	5.7
FoodTreatment	Deprived	320	53.4
	Satiated	279	46.6
SexParent	Female	245	40.9
	Male	354	59.1

owls |> modelsummary::datasummary_skim(by = c("SexParent", "FoodTreatment"))

	SexParent	FoodTreatment	Unique	Missing Pct.	Mean	SD	Min	Median	Max	Histogram
ArrivalTime	Male	Deprived	160	0	24.9	2.0	22.0	24.4	29.2
	Male	Satiated	128	0	24.7	1.8	21.8	24.5	29.1
	Female	Satiated	105	0	24.6	1.9	21.7	24.4	29.2
	Female	Deprived	107	0	24.8	2.0	22.0	24.2	29.2
SiblingNegotiation	Male	Deprived	25	0	8.7	6.5	0.0	7.0	31.0
	Male	Satiated	23	0	5.3	6.6	0.0	3.0	32.0
	Female	Satiated	22	0	4.8	6.5	0.0	1.5	26.0
	Female	Deprived	21	0	7.3	6.3	0.0	6.0	28.0
BroodSize	Male	Deprived	7	0	4.5	1.1	1.0	4.0	7.0
	Male	Satiated	7	0	4.4	1.1	1.0	5.0	7.0
	Female	Satiated	7	0	4.5	1.3	1.0	5.0	7.0
	Female	Deprived	6	0	4.0	1.1	2.0	4.0	7.0
NegPerChick	Male	Deprived	52	0	2.0	1.5	0.0	1.5	8.5
	Male	Satiated	45	0	1.2	1.5	0.0	0.6	6.7
	Female	Satiated	39	0	1.0	1.5	0.0	0.3	8.5
	Female	Deprived	44	0	1.9	1.8	0.0	1.5	8.0
		N	%
Nest	AutavauxTV	28	4.7
	Bochet	23	3.8
	Champmartin	30	5.0
	ChEsard	20	3.3
	Chevroux	10	1.7
	CorcellesFavres	12	2.0
	Etrabloz	34	5.7
	Forel	4	0.7
	Franex	26	4.3
	GDLV	10	1.7
	Gletterens	15	2.5
	Henniez	13	2.2
	Jeuss	19	3.2
	LesPlanches	17	2.8
	Lucens	29	4.8
	Lully	17	2.8
	Marnand	27	4.5
	Montet	41	6.8
	Murist	24	4.0
	Oleyes	52	8.7
	Payerne	25	4.2
	Rueyes	17	2.8
	Seiry	26	4.3
	Sevaz	4	0.7
	StAubin	23	3.8
	Trey	19	3.2
	Yvonnand	34	5.7
FoodTreatment	Deprived	320	53.4
	Satiated	279	46.6
SexParent	Female	245	40.9
	Male	354	59.1

4 Data preparation

Let start by declaring the categorical variables and random effect as factors.

## Amount of Sibling negotiation (vocalizations when parents are absent)
## Foot treatment (deprived or satiated
## Sex of parent
## Arrival time of parent
## Nest as random
## Brood size offset
owls <- owls |> mutate(
  Nest = factor(Nest),
  FoodTreatment = factor(FoodTreatment),
  SexParent = factor(SexParent),
  NCalls = SiblingNegotiation
)

5 Exploratory data analysis

As the response represents counts (the number of calls), it would make sense to start by considering a Poisson model. We could attempt to model the response as the number of calls divided by the brood size, but this would result in a response that has no natural distribution.

Instead, if we include brood size as an offset, it will standardise the effects according to brood size (similar to having divided the response by brood size), yet retain the Poisson nature of the response.

The effects of offsets, unlike regular covariates, are not estimated. Rather they are assumed to be 1 (on the link scale). This means that since Poisson uses a log link, then the offset should be of a logged version of the brood size.

Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of food treatment, sex of parent, arrival time (and various interactions) on the number of sibling negotiations. Brood size was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual nests.

Violin plots

Perhaps we could start off by exploring the main fixed effects. To mimic the log-link, we will use a log-transformed y axis. Since there may well be zeros (no calls detected), we will use a pseudo log scale). We will also include the raw data (jittered and dodged)

ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
  geom_violin()

ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
  geom_violin() +
  geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.9))

ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
  geom_violin() +
  geom_point(position = position_jitterdodge(jitter.height = 0, dodge.width = 1)) +
  scale_y_continuous(trans = scales::pseudo_log_trans())

Facetted plots

Now, a similar plot separated for each nest.

ggplot(data = owls) +
  geom_point(aes(y = NCalls, x = FoodTreatment, color = SexParent), position = position_dodge(0.5)) +
  facet_wrap(~Nest)

Smoothers

It might also be worth establishing that there is a linear relationship between the number of calls and brood size. Again, we will mimic the use of the log-link by transforming the axes.

ggplot(data = owls, aes(y = NCalls, x = BroodSize, color = SexParent)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_grid(~FoodTreatment) +
  scale_y_continuous(trans = scales::pseudo_log_trans()) +
  scale_x_log10()

`geom_smooth()` using formula = 'y ~ x'

6 Fit the model

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.

owls.rstanP <- stan_glmer(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) + (1 | Nest),
  data = owls,
  family = poisson(link = "log"),
  refresh = 0,
  iter = 5000,
  warmup = 2000,
  thin = 10,
  chains = 3,
  cores = 3
)

owls.rstanP |> prior_summary()

Priors for model 'owls.rstanP' 
------
Intercept (after predictors centered)
 ~ normal(location = 0, scale = 2.5)

Coefficients
  Specified prior:
    ~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
  Adjusted prior:
    ~ normal(location = [0,0,0], scale = [5.01,5.08,5.68])

Covariance
 ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details

This tells us:

• for the intercept (when the family is Gaussian), a normal prior with a mean of 0 and a standard deviation of 2.5 is used. The 2.5 is used for all intercepts. This is then scaled to the scale of the data by multiplying by the standard deviation of the response. When this results in values less than 2.5, it is replaced with 2.5.

2.5 / sd(owls$NCalls)

[1] 0.3747667

• for the coefficients (in this case, just the difference between strong and weak inoculation), the default prior is a normal prior centred around 0 with a standard deviation of 2.5. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the standard deviation of the numerical dummy variables for the predictor (then rounded).

model.matrix(~ FoodTreatment * SexParent + offset(log(BroodSize)), data = owls) |>
  apply(2, sd) * 1 / 2.5

                        (Intercept)               FoodTreatmentSatiated 
                          0.0000000                           0.1996977 
                      SexParentMale FoodTreatmentSatiated:SexParentMale 
                          0.1968252                           0.1760585 

Assessing priors

Lets now run with priors only so that we can explore the range of values they allow in the posteriors.

owls.rstanarmP1 <- update(owls.rstanP, prior_PD = TRUE)

owls.rstanarmP1 |>
  ggpredict(~ FoodTreatment * SexParent) |>
  plot(show_data = TRUE, jitter = c(0.25, 0)) +
  scale_y_continuous("", trans = scales::pseudo_log_trans())

You are calculating adjusted predictions on the population-level (i.e.
  `type = "fixed"`) for a *generalized* linear mixed model.
  This may produce biased estimates due to Jensen's inequality. Consider
  setting `bias_correction = TRUE` to correct for this bias.
  See also the documentation of the `bias_correction` argument.

Model uses a transformed offset term. Predictions may not be correct.
  It is recommended to fix the offset term using the `condition` argument,
  e.g. `condition = c(BroodSize = 1)`.
  You could also transform the offset variable before fitting the model
  and use `offset(BroodSize)` in the model formula.

Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Conclusions:

• we see that the range of predictions is fairly wide and the predicted means could range from a small value to a relatively large value.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

• \(\beta_0\): normal centred at 1.5 with a standard deviation of 1.5
  • mean of 1.5: since mean(log(owls$NCalls+1)) or mean(asinh(owls$NCalls/(2*1))/log(exp(1)))
  • sd of 1.5: since sd(log(owls$NCalls+1)) or sd(asinh(owls$NCalls/(2*1))/log(exp(1)))
• \(\beta_1\): normal centred at 0 with a standard deviation of 2.2
  • sd of 2.2: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\beta_2\): normal centred at 0 with a standard deviation of 2.2
  • sd of 2.2: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\beta_3\): normal centred at 0 with a standard deviation of 2.5
  • sd of 2.5: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\sigma\): exponential with rate of 1
• \(\Sigma\): decov with:
  • regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

I will also overlay the raw data for comparison.

owls.rstanarmP2 <- stan_glmer(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) + (1 | Nest),
  data = owls,
  family = poisson(link = "log"),
  prior_intercept = normal(0, 1.5, autoscale = FALSE),
  prior = normal(0, c(2.2, 2.2, 2.5), autoscale = FALSE),
  prior_aux = exponential(1),
  prior_covariance = decov(1, 1, 1, 1),
  refresh = 0,
  iter = 5000,
  prior_PD = TRUE,
  warmup = 2000,
  thin = 10,
  chains = 3,
  cores = 3
)

owls.rstanarmP2 |>
  ggpredict(~ FoodTreatment * SexParent) |>
  plot(show_data = TRUE, jitter = c(0.25, 0)) +
  scale_y_continuous("", trans = scales::pseudo_log_trans())

Model uses a transformed offset term. Predictions may not be correct.
  It is recommended to fix the offset term using the `condition` argument,
  e.g. `condition = c(BroodSize = 1)`.
  You could also transform the offset variable before fitting the model
  and use `offset(BroodSize)` in the model formula.

Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Now lets refit, conditioning on the data.

owls.rstanarmP3 <- update(owls.rstanarmP2, prior_PD = FALSE)

Plotting prior and posterior

posterior_vs_prior(owls.rstanarmP3,
  color_by = "vs", group_by = TRUE,
  facet_args = list(scales = "free_y")
)

Drawing from prior...

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

owls.form <- bf(NCalls ~ FoodTreatment*SexParent +
                    offset(log(BroodSize)) + (1|Nest),
                family=poisson(link='log'))
options(width=150)
owls.form |> get_prior(data = owls)

                               prior     class                                coef group resp dpar nlpar lb ub       source
                              (flat)         b                                                                      default
                              (flat)         b               FoodTreatmentSatiated                             (vectorized)
                              (flat)         b FoodTreatmentSatiated:SexParentMale                             (vectorized)
                              (flat)         b                       SexParentMale                             (vectorized)
 student_t(3, 0.16097150412244, 2.5) Intercept                                                                      default
                student_t(3, 0, 2.5)        sd                                                            0         default
                student_t(3, 0, 2.5)        sd                                      Nest                  0    (vectorized)
                student_t(3, 0, 2.5)        sd                           Intercept  Nest                  0    (vectorized)

options(width=80)

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

• \(\beta_0\): normal centred at 1.5 with a standard deviation of 1.5
  • mean of 1.5: since mean(log(owls$NCalls+1)) or mean(asinh(owls$NCalls/(2*1))/log(exp(1)))
  • sd of 1.5: since sd(log(owls$NCalls+1)) or sd(asinh(owls$NCalls/(2*1))/log(exp(1)))
• \(\beta_1\): normal centred at 0 with a standard deviation of 2.2
  • sd of 2.2: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\beta_2\): normal centred at 0 with a standard deviation of 2.2
  • sd of 2.2: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\beta_3\): normal centred at 0 with a standard deviation of 2.5
  • sd of 2.5: since sd(log(owls$NCalls+1))/model.matrix(~FoodTreatment*SexParent+offset(log(BroodSize)), data = owls) |> apply(2, sd)
• \(\sigma_j\): half-cauchy with parameters 0 and 5.
• \(\Sigma\): decov with:
  • regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

• half-t: as the degrees of freedom approach infinity, this will approach a half-normal
• half-cauchy: this is essentially a half-t with 1 degree of freedom
• exponential

I will also overlay the raw data for comparison.

owls |>
  group_by(FoodTreatment, SexParent) |>
  summarise(
    Mean = log(median(NCalls / BroodSize)),
    SD = log(sd(NCalls / BroodSize)),
    MAD = log(mad(NCalls / BroodSize))
  )

`summarise()` has grouped output by 'FoodTreatment'. You can override using the
`.groups` argument.

# A tibble: 4 × 5
# Groups:   FoodTreatment [2]
  FoodTreatment SexParent   Mean    SD    MAD
  <fct>         <fct>      <dbl> <dbl>  <dbl>
1 Deprived      Female     0.405 0.562  0.656
2 Deprived      Male       0.405 0.416  0.489
3 Satiated      Female    -1.23  0.399 -0.838
4 Satiated      Male      -0.511 0.431 -0.117

standist::visualize("normal(0.4, 0.5)", xlim = c(0, 20))

standist::visualize("student_t(3, 0, 2.5)",
  "cauchy(0,1)",
  xlim = c(-10, 25)
)

priors <- prior(normal(0.4, 0.7), class = "Intercept") +
  ## prior(normal(0, 2), class = 'b', coef = 'FoodTreatmentSatiated') +
  ## prior(normal(0, 2.2), class = 'b', coef = 'SexParentMale') +
  prior(normal(0, 2), class = "b") +
  prior(student_t(3, 0, 0.7), class = "sd")
## prior(cauchy(0,0.7), class = 'sd')

## priors <- prior(normal(1.8, 2), class = 'Intercept') +
##     prior(normal(0, 2.2), class = 'b', coef = 'FoodTreatmentSatiated') +
##     prior(normal(0, 2.2), class = 'b', coef = 'SexParentMale') +
##     prior(normal(0, 1), class = 'b') +
##     prior(cauchy(0,2), class = 'sd')
owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) + (1 | Nest),
  family = poisson(link = "log")
)
owls.brm2 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 0,
  seed = 123,
  control = list(adapt_delta = 0.99),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

owls.brm2 |>
  ggpredict(~ FoodTreatment * SexParent) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

owls.brm2 |>
  conditional_effects("FoodTreatment:SexParent") |>
  plot(points = TRUE)

owls.brm2 |>
  conditional_effects("FoodTreatment:SexParent") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.

owls.brm3 <- update(owls.brm2,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan",
  refresh = 0
)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

owls.brm3 |>
  ggpredict(~ FoodTreatment * SexParent) |>
  plot(show_data = TRUE)

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

owls.brm3 |>
  conditional_effects("FoodTreatment:SexParent") |>
  plot(points = TRUE) |>
  _[[1]] +
  scale_y_continuous(trans = scales::pseudo_log_trans())

Plotting prior and posterior

owls.brm3 |> get_variables()

 [1] "b_Intercept"                          
 [2] "b_FoodTreatmentSatiated"              
 [3] "b_SexParentMale"                      
 [4] "b_FoodTreatmentSatiated:SexParentMale"
 [5] "sd_Nest__Intercept"                   
 [6] "Intercept"                            
 [7] "r_Nest[AutavauxTV,Intercept]"         
 [8] "r_Nest[Bochet,Intercept]"             
 [9] "r_Nest[Champmartin,Intercept]"        
[10] "r_Nest[ChEsard,Intercept]"            
[11] "r_Nest[Chevroux,Intercept]"           
[12] "r_Nest[CorcellesFavres,Intercept]"    
[13] "r_Nest[Etrabloz,Intercept]"           
[14] "r_Nest[Forel,Intercept]"              
[15] "r_Nest[Franex,Intercept]"             
[16] "r_Nest[GDLV,Intercept]"               
[17] "r_Nest[Gletterens,Intercept]"         
[18] "r_Nest[Henniez,Intercept]"            
[19] "r_Nest[Jeuss,Intercept]"              
[20] "r_Nest[LesPlanches,Intercept]"        
[21] "r_Nest[Lucens,Intercept]"             
[22] "r_Nest[Lully,Intercept]"              
[23] "r_Nest[Marnand,Intercept]"            
[24] "r_Nest[Montet,Intercept]"             
[25] "r_Nest[Murist,Intercept]"             
[26] "r_Nest[Oleyes,Intercept]"             
[27] "r_Nest[Payerne,Intercept]"            
[28] "r_Nest[Rueyes,Intercept]"             
[29] "r_Nest[Seiry,Intercept]"              
[30] "r_Nest[Sevaz,Intercept]"              
[31] "r_Nest[StAubin,Intercept]"            
[32] "r_Nest[Trey,Intercept]"               
[33] "r_Nest[Yvonnand,Intercept]"           
[34] "prior_Intercept"                      
[35] "prior_b"                              
[36] "prior_sd_Nest"                        
[37] "lprior"                               
[38] "lp__"                                 
[39] "accept_stat__"                        
[40] "stepsize__"                           
[41] "treedepth__"                          
[42] "n_leapfrog__"                         
[43] "divergent__"                          
[44] "energy__"                             

owls.brm3 |>
  hypothesis("FoodTreatmentSatiated=0") |>
  plot()

owls.brm3 |>
  hypothesis("SexParentMale=0") |>
  plot()

owls.brm3 |> SUYR_prior_and_posterior()

owls.brm3 |>
  posterior_samples() |>
  dplyr::select(-`lp__`) |>
  pivot_longer(everything(), names_to = "key") |>
  filter(!str_detect(key, "^r")) |>
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    ## Class = ifelse(str_detect(key, 'Intercept'),  'Intercept',
    ##         ifelse(str_detect(key, 'b'),  'b', 'sigma')),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_FoodTreatment.*|b_SexParent.*|prior_b") ~ "TREATMENT",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) |>
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
  facet_wrap(~Class, scales = "free")

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

Random intercept/slope model

While we are here, we might like to explore a random intercept/slope model

priors <- prior(normal(0.4, 0.7), class = "Intercept") +
  ## prior(normal(0, 2.2), class = 'b', coef = 'FoodTreatmentSatiated') +
  ## prior(normal(0, 2.2), class = 'b', coef = 'SexParentMale') +
  prior(normal(0, 2), class = "b") +
  prior(cauchy(0, 0.7), class = "sd") +
  prior(cauchy(0, 0.7), class = "sd", coef = "Intercept", group = "Nest") +
  prior(lkj_corr_cholesky(1), class = "cor")

## 130 seconds

owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) +
    (FoodTreatment * SexParent | Nest),
  family = poisson(link = "log")
)
owls.brm4 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "yes",
  ## sample_prior = 'only',
  iter = 5000,
  warmup = 2500,
  chains = 3,
  cores = 3,
  thin = 10,
  refresh = 100,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

owls.brm4 |> get_variables()

  [1] "b_Intercept"                                                         
  [2] "b_FoodTreatmentSatiated"                                             
  [3] "b_SexParentMale"                                                     
  [4] "b_FoodTreatmentSatiated:SexParentMale"                               
  [5] "sd_Nest__Intercept"                                                  
  [6] "sd_Nest__FoodTreatmentSatiated"                                      
  [7] "sd_Nest__SexParentMale"                                              
  [8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"                        
  [9] "cor_Nest__Intercept__FoodTreatmentSatiated"                          
 [10] "cor_Nest__Intercept__SexParentMale"                                  
 [11] "cor_Nest__FoodTreatmentSatiated__SexParentMale"                      
 [12] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"            
 [13] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
 [14] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"        
 [15] "Intercept"                                                           
 [16] "r_Nest[AutavauxTV,Intercept]"                                        
 [17] "r_Nest[Bochet,Intercept]"                                            
 [18] "r_Nest[Champmartin,Intercept]"                                       
 [19] "r_Nest[ChEsard,Intercept]"                                           
 [20] "r_Nest[Chevroux,Intercept]"                                          
 [21] "r_Nest[CorcellesFavres,Intercept]"                                   
 [22] "r_Nest[Etrabloz,Intercept]"                                          
 [23] "r_Nest[Forel,Intercept]"                                             
 [24] "r_Nest[Franex,Intercept]"                                            
 [25] "r_Nest[GDLV,Intercept]"                                              
 [26] "r_Nest[Gletterens,Intercept]"                                        
 [27] "r_Nest[Henniez,Intercept]"                                           
 [28] "r_Nest[Jeuss,Intercept]"                                             
 [29] "r_Nest[LesPlanches,Intercept]"                                       
 [30] "r_Nest[Lucens,Intercept]"                                            
 [31] "r_Nest[Lully,Intercept]"                                             
 [32] "r_Nest[Marnand,Intercept]"                                           
 [33] "r_Nest[Montet,Intercept]"                                            
 [34] "r_Nest[Murist,Intercept]"                                            
 [35] "r_Nest[Oleyes,Intercept]"                                            
 [36] "r_Nest[Payerne,Intercept]"                                           
 [37] "r_Nest[Rueyes,Intercept]"                                            
 [38] "r_Nest[Seiry,Intercept]"                                             
 [39] "r_Nest[Sevaz,Intercept]"                                             
 [40] "r_Nest[StAubin,Intercept]"                                           
 [41] "r_Nest[Trey,Intercept]"                                              
 [42] "r_Nest[Yvonnand,Intercept]"                                          
 [43] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"                            
 [44] "r_Nest[Bochet,FoodTreatmentSatiated]"                                
 [45] "r_Nest[Champmartin,FoodTreatmentSatiated]"                           
 [46] "r_Nest[ChEsard,FoodTreatmentSatiated]"                               
 [47] "r_Nest[Chevroux,FoodTreatmentSatiated]"                              
 [48] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"                       
 [49] "r_Nest[Etrabloz,FoodTreatmentSatiated]"                              
 [50] "r_Nest[Forel,FoodTreatmentSatiated]"                                 
 [51] "r_Nest[Franex,FoodTreatmentSatiated]"                                
 [52] "r_Nest[GDLV,FoodTreatmentSatiated]"                                  
 [53] "r_Nest[Gletterens,FoodTreatmentSatiated]"                            
 [54] "r_Nest[Henniez,FoodTreatmentSatiated]"                               
 [55] "r_Nest[Jeuss,FoodTreatmentSatiated]"                                 
 [56] "r_Nest[LesPlanches,FoodTreatmentSatiated]"                           
 [57] "r_Nest[Lucens,FoodTreatmentSatiated]"                                
 [58] "r_Nest[Lully,FoodTreatmentSatiated]"                                 
 [59] "r_Nest[Marnand,FoodTreatmentSatiated]"                               
 [60] "r_Nest[Montet,FoodTreatmentSatiated]"                                
 [61] "r_Nest[Murist,FoodTreatmentSatiated]"                                
 [62] "r_Nest[Oleyes,FoodTreatmentSatiated]"                                
 [63] "r_Nest[Payerne,FoodTreatmentSatiated]"                               
 [64] "r_Nest[Rueyes,FoodTreatmentSatiated]"                                
 [65] "r_Nest[Seiry,FoodTreatmentSatiated]"                                 
 [66] "r_Nest[Sevaz,FoodTreatmentSatiated]"                                 
 [67] "r_Nest[StAubin,FoodTreatmentSatiated]"                               
 [68] "r_Nest[Trey,FoodTreatmentSatiated]"                                  
 [69] "r_Nest[Yvonnand,FoodTreatmentSatiated]"                              
 [70] "r_Nest[AutavauxTV,SexParentMale]"                                    
 [71] "r_Nest[Bochet,SexParentMale]"                                        
 [72] "r_Nest[Champmartin,SexParentMale]"                                   
 [73] "r_Nest[ChEsard,SexParentMale]"                                       
 [74] "r_Nest[Chevroux,SexParentMale]"                                      
 [75] "r_Nest[CorcellesFavres,SexParentMale]"                               
 [76] "r_Nest[Etrabloz,SexParentMale]"                                      
 [77] "r_Nest[Forel,SexParentMale]"                                         
 [78] "r_Nest[Franex,SexParentMale]"                                        
 [79] "r_Nest[GDLV,SexParentMale]"                                          
 [80] "r_Nest[Gletterens,SexParentMale]"                                    
 [81] "r_Nest[Henniez,SexParentMale]"                                       
 [82] "r_Nest[Jeuss,SexParentMale]"                                         
 [83] "r_Nest[LesPlanches,SexParentMale]"                                   
 [84] "r_Nest[Lucens,SexParentMale]"                                        
 [85] "r_Nest[Lully,SexParentMale]"                                         
 [86] "r_Nest[Marnand,SexParentMale]"                                       
 [87] "r_Nest[Montet,SexParentMale]"                                        
 [88] "r_Nest[Murist,SexParentMale]"                                        
 [89] "r_Nest[Oleyes,SexParentMale]"                                        
 [90] "r_Nest[Payerne,SexParentMale]"                                       
 [91] "r_Nest[Rueyes,SexParentMale]"                                        
 [92] "r_Nest[Seiry,SexParentMale]"                                         
 [93] "r_Nest[Sevaz,SexParentMale]"                                         
 [94] "r_Nest[StAubin,SexParentMale]"                                       
 [95] "r_Nest[Trey,SexParentMale]"                                          
 [96] "r_Nest[Yvonnand,SexParentMale]"                                      
 [97] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"              
 [98] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"                  
 [99] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"             
[100] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"                 
[101] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"                
[102] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"         
[103] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"                
[104] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"                   
[105] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"                  
[106] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"                    
[107] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"              
[108] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"                 
[109] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"                   
[110] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"             
[111] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"                  
[112] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"                   
[113] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"                 
[114] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"                  
[115] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"                  
[116] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"                  
[117] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"                 
[118] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"                  
[119] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"                   
[120] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"                   
[121] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"                 
[122] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"                    
[123] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"                
[124] "prior_Intercept"                                                     
[125] "prior_b"                                                             
[126] "prior_sd_Nest__Intercept"                                            
[127] "prior_sd_Nest__FoodTreatmentSatiated"                                
[128] "prior_sd_Nest__SexParentMale"                                        
[129] "prior_sd_Nest__FoodTreatmentSatiated:SexParentMale"                  
[130] "prior_cor_Nest"                                                      
[131] "lprior"                                                              
[132] "lp__"                                                                
[133] "accept_stat__"                                                       
[134] "stepsize__"                                                          
[135] "treedepth__"                                                         
[136] "n_leapfrog__"                                                        
[137] "divergent__"                                                         
[138] "energy__"                                                            

owls.brm4 |>
  hypothesis("FoodTreatmentSatiated=0") |>
  plot()

owls.brm4 |>
  hypothesis("SexParentMale=0") |>
  plot()

owls.brm4 |> SUYR_prior_and_posterior()

owls.brm4 |>
  posterior_samples() |>
  dplyr::select(-`lp__`) |>
  pivot_longer(everything(), names_to = "key") |>
  filter(!str_detect(key, "^r")) |>
  mutate(
    Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
    ## Class = ifelse(str_detect(key, 'Intercept'),  'Intercept',
    ##         ifelse(str_detect(key, 'b'),  'b', 'sigma')),
    Class = case_when(
      str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
      str_detect(key, "b_FoodTreatment.*|prior_b_FoodTreatment.*") &
        !str_detect(key, ".*:.*") ~ "FoodTreatment",
      str_detect(key, "b_SexParent.*|prior_b_SexParent.*") &
        !str_detect(key, ".*\\:.*") ~ "SexParent",
      str_detect(key, ".*\\:.*|prior_b_.*\\:.*") ~ "Interaction",
      str_detect(key, "sd") ~ "sd",
      str_detect(key, "^cor|prior_cor") ~ "cor",
      str_detect(key, "sigma") ~ "sigma"
    ),
    Par = str_replace(key, "b_", "")
  ) |>
  ggplot(aes(x = Type, y = value, color = Par)) +
  stat_pointinterval(position = position_dodge()) +
  facet_wrap(~Class, scales = "free")

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

Compare models

We can compare the two models using LOO

(l.1 <- owls.brm3 |> loo::loo())

Warning: Found 3 observations with a pareto_k > 0.7 in model 'owls.brm3'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.

Computed from 750 by 599 log-likelihood matrix.

         Estimate    SE
elpd_loo  -2636.0  85.8
p_loo       145.0  11.2
looic      5271.9 171.6
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.3]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.65]   (good)     596   99.5%   49      
   (0.65, 1]   (bad)        3    0.5%   <NA>    
    (1, Inf)   (very bad)   0    0.0%   <NA>    
See help('pareto-k-diagnostic') for details.

(l.2 <- owls.brm4 |> loo::loo())

Warning: Found 43 observations with a pareto_k > 0.65 in model 'owls.brm4'. We
recommend to run more iterations to get at least about 2200 posterior draws to
improve LOO-CV approximation accuracy.

Computed from 750 by 599 log-likelihood matrix.

         Estimate    SE
elpd_loo  -2427.1  74.5
p_loo       296.4  19.6
looic      4854.2 149.1
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.65]   (good)     556   92.8%   32      
   (0.65, 1]   (bad)       35    5.8%   <NA>    
    (1, Inf)   (very bad)   8    1.3%   <NA>    
See help('pareto-k-diagnostic') for details.

loo::loo_compare(l.1, l.2)

          elpd_diff se_diff
owls.brm4    0.0       0.0 
owls.brm3 -208.9      47.1 

Although there is not substantially more support for the more complex random intercept/slope model over the simpler random intercept only model, we might go with the more complex model anyway.

7 MCMC sampling diagnostics

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

pars <- owls.brm4 |> get_variables()
pars <- pars |>
  str_extract("^b.Intercept|^b_FoodTreatment.*|^b_SexParent.*|[sS]igma|^sd.*") |>
  na.omit()
pars

[1] "b_Intercept"                                 
[2] "b_FoodTreatmentSatiated"                     
[3] "b_SexParentMale"                             
[4] "b_FoodTreatmentSatiated:SexParentMale"       
[5] "sd_Nest__Intercept"                          
[6] "sd_Nest__FoodTreatmentSatiated"              
[7] "sd_Nest__SexParentMale"                      
[8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"
attr(,"na.action")
  [1]   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26
 [19]  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44
 [37]  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62
 [55]  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80
 [73]  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98
 [91]  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
[109] 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
[127] 135 136 137 138
attr(,"class")
[1] "omit"

owls.brm4 |> mcmc_plot(type = "trace", variable = pars)

No divergences to plot.

## OR
owls.brm4 |> mcmc_plot(
  type = "trace",
  variable = "^b.Intercept|^b_FoodTreatment.*|^b_SexParent.*|[sS]igma|^sd.*",
  regex = TRUE
)

No divergences to plot.

The chains appear well mixed and very similar

• acf_bar (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

owls.brm4 |> mcmc_plot(type = "acf_bar", variable = pars)

## OR
owls.brm4 |> mcmc_plot(
  type = "acf_bar",
  variable = "^b.Intercept|^b_FoodTreatment.*|^b_SexParent.*|[sS]igma|^sd.*",
  regex = TRUE
)

There is no evidence of auto-correlation in the MCMC samples

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

owls.brm4 |> mcmc_plot(type = "rhat_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

owls.brm4 |> mcmc_plot(type = "neff_hist")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics

owls.brm4 |> mcmc_plot(type = "combo", pars = pars)

Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

owls.brm4 |> mcmc_plot(type = "violin", pars = pars)

Warning: Argument 'pars' is deprecated. Please use 'variable' instead.

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

owls.brm4 |> get_variables()

  [1] "b_Intercept"                                                         
  [2] "b_FoodTreatmentSatiated"                                             
  [3] "b_SexParentMale"                                                     
  [4] "b_FoodTreatmentSatiated:SexParentMale"                               
  [5] "sd_Nest__Intercept"                                                  
  [6] "sd_Nest__FoodTreatmentSatiated"                                      
  [7] "sd_Nest__SexParentMale"                                              
  [8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"                        
  [9] "cor_Nest__Intercept__FoodTreatmentSatiated"                          
 [10] "cor_Nest__Intercept__SexParentMale"                                  
 [11] "cor_Nest__FoodTreatmentSatiated__SexParentMale"                      
 [12] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"            
 [13] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
 [14] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"        
 [15] "Intercept"                                                           
 [16] "r_Nest[AutavauxTV,Intercept]"                                        
 [17] "r_Nest[Bochet,Intercept]"                                            
 [18] "r_Nest[Champmartin,Intercept]"                                       
 [19] "r_Nest[ChEsard,Intercept]"                                           
 [20] "r_Nest[Chevroux,Intercept]"                                          
 [21] "r_Nest[CorcellesFavres,Intercept]"                                   
 [22] "r_Nest[Etrabloz,Intercept]"                                          
 [23] "r_Nest[Forel,Intercept]"                                             
 [24] "r_Nest[Franex,Intercept]"                                            
 [25] "r_Nest[GDLV,Intercept]"                                              
 [26] "r_Nest[Gletterens,Intercept]"                                        
 [27] "r_Nest[Henniez,Intercept]"                                           
 [28] "r_Nest[Jeuss,Intercept]"                                             
 [29] "r_Nest[LesPlanches,Intercept]"                                       
 [30] "r_Nest[Lucens,Intercept]"                                            
 [31] "r_Nest[Lully,Intercept]"                                             
 [32] "r_Nest[Marnand,Intercept]"                                           
 [33] "r_Nest[Montet,Intercept]"                                            
 [34] "r_Nest[Murist,Intercept]"                                            
 [35] "r_Nest[Oleyes,Intercept]"                                            
 [36] "r_Nest[Payerne,Intercept]"                                           
 [37] "r_Nest[Rueyes,Intercept]"                                            
 [38] "r_Nest[Seiry,Intercept]"                                             
 [39] "r_Nest[Sevaz,Intercept]"                                             
 [40] "r_Nest[StAubin,Intercept]"                                           
 [41] "r_Nest[Trey,Intercept]"                                              
 [42] "r_Nest[Yvonnand,Intercept]"                                          
 [43] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"                            
 [44] "r_Nest[Bochet,FoodTreatmentSatiated]"                                
 [45] "r_Nest[Champmartin,FoodTreatmentSatiated]"                           
 [46] "r_Nest[ChEsard,FoodTreatmentSatiated]"                               
 [47] "r_Nest[Chevroux,FoodTreatmentSatiated]"                              
 [48] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"                       
 [49] "r_Nest[Etrabloz,FoodTreatmentSatiated]"                              
 [50] "r_Nest[Forel,FoodTreatmentSatiated]"                                 
 [51] "r_Nest[Franex,FoodTreatmentSatiated]"                                
 [52] "r_Nest[GDLV,FoodTreatmentSatiated]"                                  
 [53] "r_Nest[Gletterens,FoodTreatmentSatiated]"                            
 [54] "r_Nest[Henniez,FoodTreatmentSatiated]"                               
 [55] "r_Nest[Jeuss,FoodTreatmentSatiated]"                                 
 [56] "r_Nest[LesPlanches,FoodTreatmentSatiated]"                           
 [57] "r_Nest[Lucens,FoodTreatmentSatiated]"                                
 [58] "r_Nest[Lully,FoodTreatmentSatiated]"                                 
 [59] "r_Nest[Marnand,FoodTreatmentSatiated]"                               
 [60] "r_Nest[Montet,FoodTreatmentSatiated]"                                
 [61] "r_Nest[Murist,FoodTreatmentSatiated]"                                
 [62] "r_Nest[Oleyes,FoodTreatmentSatiated]"                                
 [63] "r_Nest[Payerne,FoodTreatmentSatiated]"                               
 [64] "r_Nest[Rueyes,FoodTreatmentSatiated]"                                
 [65] "r_Nest[Seiry,FoodTreatmentSatiated]"                                 
 [66] "r_Nest[Sevaz,FoodTreatmentSatiated]"                                 
 [67] "r_Nest[StAubin,FoodTreatmentSatiated]"                               
 [68] "r_Nest[Trey,FoodTreatmentSatiated]"                                  
 [69] "r_Nest[Yvonnand,FoodTreatmentSatiated]"                              
 [70] "r_Nest[AutavauxTV,SexParentMale]"                                    
 [71] "r_Nest[Bochet,SexParentMale]"                                        
 [72] "r_Nest[Champmartin,SexParentMale]"                                   
 [73] "r_Nest[ChEsard,SexParentMale]"                                       
 [74] "r_Nest[Chevroux,SexParentMale]"                                      
 [75] "r_Nest[CorcellesFavres,SexParentMale]"                               
 [76] "r_Nest[Etrabloz,SexParentMale]"                                      
 [77] "r_Nest[Forel,SexParentMale]"                                         
 [78] "r_Nest[Franex,SexParentMale]"                                        
 [79] "r_Nest[GDLV,SexParentMale]"                                          
 [80] "r_Nest[Gletterens,SexParentMale]"                                    
 [81] "r_Nest[Henniez,SexParentMale]"                                       
 [82] "r_Nest[Jeuss,SexParentMale]"                                         
 [83] "r_Nest[LesPlanches,SexParentMale]"                                   
 [84] "r_Nest[Lucens,SexParentMale]"                                        
 [85] "r_Nest[Lully,SexParentMale]"                                         
 [86] "r_Nest[Marnand,SexParentMale]"                                       
 [87] "r_Nest[Montet,SexParentMale]"                                        
 [88] "r_Nest[Murist,SexParentMale]"                                        
 [89] "r_Nest[Oleyes,SexParentMale]"                                        
 [90] "r_Nest[Payerne,SexParentMale]"                                       
 [91] "r_Nest[Rueyes,SexParentMale]"                                        
 [92] "r_Nest[Seiry,SexParentMale]"                                         
 [93] "r_Nest[Sevaz,SexParentMale]"                                         
 [94] "r_Nest[StAubin,SexParentMale]"                                       
 [95] "r_Nest[Trey,SexParentMale]"                                          
 [96] "r_Nest[Yvonnand,SexParentMale]"                                      
 [97] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"              
 [98] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"                  
 [99] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"             
[100] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"                 
[101] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"                
[102] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"         
[103] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"                
[104] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"                   
[105] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"                  
[106] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"                    
[107] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"              
[108] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"                 
[109] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"                   
[110] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"             
[111] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"                  
[112] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"                   
[113] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"                 
[114] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"                  
[115] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"                  
[116] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"                  
[117] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"                 
[118] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"                  
[119] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"                   
[120] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"                   
[121] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"                 
[122] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"                    
[123] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"                
[124] "prior_Intercept"                                                     
[125] "prior_b"                                                             
[126] "prior_sd_Nest__Intercept"                                            
[127] "prior_sd_Nest__FoodTreatmentSatiated"                                
[128] "prior_sd_Nest__SexParentMale"                                        
[129] "prior_sd_Nest__FoodTreatmentSatiated:SexParentMale"                  
[130] "prior_cor_Nest"                                                      
[131] "lprior"                                                              
[132] "lp__"                                                                
[133] "accept_stat__"                                                       
[134] "stepsize__"                                                          
[135] "treedepth__"                                                         
[136] "n_leapfrog__"                                                        
[137] "divergent__"                                                         
[138] "energy__"                                                            

pars <- owls.brm4 |> get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") |> na.omit()

owls.brm4$fit |>
  stan_trace(pars = pars)

The chains appear well mixed and very similar

• stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

owls.brm4$fit |>
  stan_ac(pars = pars)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

owls.brm4$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

owls.brm4$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

owls.brm4$fit |>
  stan_dens(separate_chains = TRUE, pars = pars)

ggmcmc

The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

## owls.ggs <- owls.brm3 |> ggs(burnin = FALSE, inc_warmup = FALSE)
## owls.ggs |> ggs_traceplot()

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

## ggs_autocorrelation(owls.ggs)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

## ggs_Rhat(owls.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

## ggs_effective(owls.ggs)

Ratios all very high.

More diagnostics

## ggs_crosscorrelation(owls.ggs)

## ggs_grb(owls.ggs)

8 Model validation

brms

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed on top of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

owls.brm4 |> pp_check(type = "dens_overlay", ndraws = 100)

The model draws appear to differ substantially from the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

owls.brm4 |> pp_check(type = "error_scatter_avg")

Using all posterior draws for ppc type 'error_scatter_avg' by default.

This is not really interpretable

• intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

owls.brm4 |> pp_check(group = "Nest", type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

Warning: The following arguments were unrecognized and ignored: group

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

# library(shinystan)
# launch_shinystan(owls.brm2)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

## preds <- owls.brm4 |> posterior_predict(nsamples = 250,  summary = FALSE)
## owls.resids <- createDHARMa(simulatedResponse = t(preds),
##                             observedResponse = owls$NCalls,
##                             fittedPredictedResponse = apply(preds, 2, median),
##                             integerResponse = TRUE)
## plot(owls.resids)

owls.resids <- make_brms_dharma_res(owls.brm4, integerResponse = TRUE)
wrap_elements(~ testUniformity(owls.resids)) +
  wrap_elements(~ plotResiduals(owls.resids, form = factor(rep(1, nrow(owls))))) +
  wrap_elements(~ plotResiduals(owls.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(owls.resids))

Conclusions:

• the model does not appear to be a very good fit
• the Q-Q plot deviates substantially from a straight line
• there are outliers

Perhaps we should specifically explore zero-inflation.

owls.resids |> testZeroInflation()

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 2.3734, p-value = 0.002667
alternative hypothesis: two.sided

Conclusions:

• there is strong evidence of zero-inflation

The data were collected at various times throughout the night. It is possible that this could lead to patterns of dependency that are not already accounted for. For example, perhaps observations that are collected at similar time of the night (within a given nest) have more similar residuals than those at very different time of the night. We can explore whether there are any temporal auto-correlation patterns.

owls.resids |> testTemporalAutocorrelation(time = owls$ArrivalTime)

Error in testTemporalAutocorrelation(owls.resids, time = owls$ArrivalTime): testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use the recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testTemporalAutocorrelation.

## owls.resid1 <- owls.resids |> recalculateResiduals(group=interaction(owls$ArrivalTime,  owls$Nest),  aggregateBy = mean)
## owls.resid1 <- owls.resids |> recalculateResiduals(group=interaction(owls$ArrivalTime,  owls$Nest),  aggregateBy = sum)
## resids1 <- owls.resids |> recalculateResiduals(group = interaction(owls$ArrivalTime), aggregateBy = sum)
resids1 <- owls.resids |> recalculateResiduals(group = interaction(owls$ArrivalTime), aggregateBy = mean)
resids1 |> testTemporalAutocorrelation(time = unique(owls$ArrivalTime))

    Durbin-Watson test

data:  simulationOutput$scaledResiduals ~ 1
DW = 2.008, p-value = 0.943
alternative hypothesis: true autocorrelation is not 0

library(geoR)
autocor_check(owls, owls.brm4,
  variable = "ArrivalTime",
  grouping = "Nest",
  n.sim = 250
)

variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin
variog: computing omnidirectional variogram
variog: co-locatted data found, adding one bin at the origin

Error in loadNamespace(x): there is no package called 'reshape'

Conclusions:

• there is no evidence of temporal auto-correlation

Conclusions:

• there is evidence that the model does not fit that well. It is evidently zero inflated and possibly also over-dispersed.
• it would seem that a zero-inflated Poisson or even a zero-inflated Negative Binomial would be a sensible next step.
• zero-inflated models cannot be fit in glmer(), so we will proceed with glmmTMB() only.

9 Model refit and validation

Fit model 1 (Poisson)

Fit

owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) +
    (FoodTreatment * SexParent | Nest),
  zi ~ 1,
  family = zero_inflated_poisson(link = "log")
)
priors <- prior(normal(0.4, 0.7), class = "Intercept") +
  ## prior(normal(0, 2.2), class = 'b', coef = 'FoodTreatmentSatiated') +
  ## prior(normal(0, 2.2), class = 'b', coef = 'SexParentMale') +
  prior(normal(0, 2), class = "b") +
  prior(student_t(3, 0, 0.7), class = "sd") +
  prior(student_t(3, 0, 0.7), class = "sd", coef = "Intercept", group = "Nest") +
  prior(lkj_corr_cholesky(1), class = "cor") +
  prior(logistic(0, 1), class = "Intercept", dpar = "zi")

## 600 seconds

owls.brm5 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "yes",
  iter = 10000,
  warmup = 5000,
  thin = 10,
  chains = 3, cores = 3,
  refresh = 0,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)

Compiling Stan program...

Start sampling

Model validation

owls.brm5 |> SUYR_prior_and_posterior()

owls.brm5 |> pp_check(type = "dens_overlay", ndraws = 100)

owls.brm5 |> pp_check(type = "dens_overlay", ndraws = 100) +
  scale_x_log10()

Warning in scale_x_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.

Warning: Removed 15040 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 156 rows containing non-finite outside the scale range
(`stat_density()`).

owls.resids <- make_brms_dharma_res(owls.brm5, integerResponse = TRUE)
wrap_elements(~ testUniformity(owls.resids)) +
  wrap_elements(~ plotResiduals(owls.resids, form = factor(rep(1, nrow(owls))))) +
  wrap_elements(~ plotResiduals(owls.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(owls.resids))

Error in graphics::par(old_gp) : 
  invalid value specified for graphical parameter "pin"

preds <- owls.brm5 |> posterior_predict(nsamples = 250, summary = FALSE)

Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.

owls.resids |> testZeroInflation()

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 0.98847, p-value = 0.968
alternative hypothesis: two.sided

owls.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.052094, p-value = 0.196
alternative hypothesis: two.sided

owls.resids |> testUniformity()

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.069154, p-value = 0.0065
alternative hypothesis: two-sided

owls.resids |> testQuantiles()

    Test for location of quantiles via qgam

data:  res
p-value = 0.01301
alternative hypothesis: both

owls.resids |> testResiduals()

$uniformity

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.069154, p-value = 0.0065
alternative hypothesis: two-sided


$dispersion

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.052094, p-value = 0.196
alternative hypothesis: two.sided


$outliers

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 2, observations = 599, p-value = 0.1905
alternative hypothesis: true probability of success is not equal to 0.001332445
95 percent confidence interval:
 0.0004046121 0.0120087614
sample estimates:
frequency of outliers (expected: 0.00133244503664224 ) 
                                           0.003338898 

Fit model 2 (Poisson)

Fit

priors <- prior(normal(1.5, 1.5), class = "Intercept") +
  prior(normal(0, 2.2), class = "b", coef = "FoodTreatmentSatiated") +
  prior(normal(0, 2.2), class = "b", coef = "SexParentMale") +
  prior(normal(0, 2.5), class = "b") +
  prior(cauchy(0, 1), class = "sd") +
  prior(lkj_corr_cholesky(1), class = "L") +
  prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
  prior(normal(0, 1), class = "b", dpar = "zi")
owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) +
    (FoodTreatment * SexParent | Nest),
  zi ~ FoodTreatment * SexParent,
  family = zero_inflated_poisson(link = "log")
)

owls.brm6 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  thin = 10,
  chains = 3,
  refresh = 0,
  backend = "rstan",
  cores = 3
)

Compiling Stan program...

Start sampling

Model validation

preds <- owls.brm6 |> posterior_predict(ndraws = 250, summary = FALSE)
owls.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = owls$NCalls,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)

DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

owls.resids |> testZeroInflation()

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 1.0143, p-value = 0.928
alternative hypothesis: two.sided

owls.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.4344, p-value = 0.032
alternative hypothesis: two.sided

owls.resids |> testUniformity()

DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.086635, p-value = 0.0002489
alternative hypothesis: two-sided

owls.resids |> testQuantiles()

    Test for location of quantiles via qgam

data:  res
p-value = 4.324e-06
alternative hypothesis: both

owls.resids |> testResiduals()

DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

$uniformity

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.086635, p-value = 0.0002489
alternative hypothesis: two-sided


$dispersion

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 1.4344, p-value = 0.032
alternative hypothesis: two.sided


$outliers

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 26, observations = 599, p-value =
7.823e-12
alternative hypothesis: true probability of success is not equal to 0.007968127
95 percent confidence interval:
 0.02854681 0.06295488
sample estimates:
frequency of outliers (expected: 0.00796812749003984 ) 
                                            0.04340568 

Fit model 3 (Negative Binomial)

Fit

priors <- prior(normal(1.5, 1.5), class = "Intercept") +
  prior(normal(0, 2.2), class = "b", coef = "FoodTreatmentSatiated") +
  prior(normal(0, 2.2), class = "b", coef = "SexParentMale") +
  prior(normal(0, 2.5), class = "b") +
  prior(cauchy(0, 1), class = "sd") +
  prior(lkj_corr_cholesky(1), class = "L") +
  prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
  ## prior(normal(0,1), class='b', dpar='zi') +
  prior(gamma(0.01, 0.01), class = "shape")
owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) +
    (FoodTreatment * SexParent | Nest),
  zi ~ 1,
  family = zero_inflated_negbinomial(link = "log")
)

owls.brm7 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "yes",
  iter = 5000,
  warmup = 2500,
  thin = 10,
  chains = 3,
  refresh = 0,
  backend = "rstan",
  cores = 3
)

Compiling Stan program...

Start sampling

Model validation

preds <- owls.brm7 |> posterior_predict(ndraws = 250, summary = FALSE)
owls.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = owls$NCalls,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)

owls.resids |> testZeroInflation()

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 1.054, p-value = 0.632
alternative hypothesis: two.sided

owls.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.69166, p-value = 0.016
alternative hypothesis: two.sided

owls.resids |> testUniformity()

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.035264, p-value = 0.4457
alternative hypothesis: two-sided

owls.resids |> testQuantiles()

    Test for location of quantiles via qgam

data:  res
p-value = 0.0008583
alternative hypothesis: both

owls.resids |> testResiduals()

$uniformity

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.035264, p-value = 0.4457
alternative hypothesis: two-sided


$dispersion

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.69166, p-value = 0.016
alternative hypothesis: two.sided


$outliers

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 2, observations = 599, p-value = 0.3488
alternative hypothesis: true probability of success is not equal to 0.007968127
95 percent confidence interval:
 0.0004046121 0.0120087614
sample estimates:
frequency of outliers (expected: 0.00796812749003984 ) 
                                           0.003338898 

Fit model 4 (Negative Binomial)

Fit

owls.form <- bf(
  NCalls ~ FoodTreatment * SexParent +
    offset(log(BroodSize)) +
    (FoodTreatment * SexParent | Nest),
  zi ~ FoodTreatment * SexParent,
  family = zero_inflated_negbinomial(link = "log")
)
priors <- prior(normal(0.4, 0.7), class = "Intercept") +
  ## prior(normal(0, 2), class = 'b', coef = 'FoodTreatmentSatiated') +
  ## prior(normal(0, 2), class = 'b', coef = 'SexParentMale') +
  prior(normal(0, 2), class = "b") +
  # prior(cauchy(0,2), class = 'sd') +
  prior(student_t(3, 0, 0.7), class = "sd") +
  prior(lkj_corr_cholesky(1), class = "L") +
  prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
  prior(normal(0, 1), class = "b", dpar = "zi") +
  prior(gamma(0.01, 0.01), class = "shape")
## 424
owls.brm8 <- brm(owls.form,
  data = owls,
  prior = priors,
  sample_prior = "yes",
  iter = 10000,
  warmup = 5000,
  thin = 10,
  chains = 3, cores = 3,
  refresh = 0,
  backend = "rstan",
  control = list(adapt_delta = 0.99)
)

Compiling Stan program...

Start sampling

Model validation

## owls.brm8 |> SUYR_prior_and_posterior()


owls.brm8 |> pp_check(type = "dens_overlay", ndraws = 100)

owls.brm8 |> pp_check(type = "dens_overlay", ndraws = 100) +
  scale_x_log10()

Warning in scale_x_log10(): log-10 transformation introduced infinite values.
log-10 transformation introduced infinite values.

Warning: Removed 15196 rows containing non-finite outside the scale range
(`stat_density()`).

Warning: Removed 156 rows containing non-finite outside the scale range
(`stat_density()`).

owls.resids <- make_brms_dharma_res(owls.brm8, integerResponse = TRUE)
wrap_elements(~ testUniformity(owls.resids)) +
  wrap_elements(~ plotResiduals(owls.resids, form = factor(rep(1, nrow(owls))))) +
  wrap_elements(~ plotResiduals(owls.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(owls.resids))

Error in graphics::par(old_gp) : 
  invalid value specified for graphical parameter "pin"

owls.resids |> testZeroInflation()

    DHARMa zero-inflation test via comparison to expected zeros with
    simulation under H0 = fitted model

data:  simulationOutput
ratioObsSim = 0.99275, p-value = 0.9827
alternative hypothesis: two.sided

owls.resids |> testDispersion()

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.070108, p-value = 0.07333
alternative hypothesis: two.sided

owls.resids |> testUniformity()

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.058502, p-value = 0.03314
alternative hypothesis: two-sided

owls.resids |> testQuantiles()

Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
: Fitting terminated with step failure - check results carefully

    Test for location of quantiles via qgam

data:  res
p-value = 0.1359
alternative hypothesis: both

owls.resids |> testResiduals()

$uniformity

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.058502, p-value = 0.03314
alternative hypothesis: two-sided


$dispersion

    DHARMa nonparametric dispersion test via sd of residuals fitted vs.
    simulated

data:  simulationOutput
dispersion = 0.070108, p-value = 0.07333
alternative hypothesis: two.sided


$outliers

    DHARMa outlier test based on exact binomial test with approximate
    expectations

data:  simulationOutput
outliers at both margin(s) = 2, observations = 599, p-value = 0.1905
alternative hypothesis: true probability of success is not equal to 0.001332445
95 percent confidence interval:
 0.0004046121 0.0120087614
sample estimates:
frequency of outliers (expected: 0.00133244503664224 ) 
                                           0.003338898 

Comparisons

(l.1 <- loo::loo(owls.brm4))

Warning: Found 43 observations with a pareto_k > 0.65 in model 'owls.brm4'. We
recommend to run more iterations to get at least about 2200 posterior draws to
improve LOO-CV approximation accuracy.

Computed from 750 by 599 log-likelihood matrix.

         Estimate    SE
elpd_loo  -2427.1  74.5
p_loo       296.4  19.6
looic      4854.2 149.1
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.65]   (good)     556   92.8%   32      
   (0.65, 1]   (bad)       35    5.8%   <NA>    
    (1, Inf)   (very bad)   8    1.3%   <NA>    
See help('pareto-k-diagnostic') for details.

(l.2 <- loo::loo(owls.brm5))

Warning: Found 26 observations with a pareto_k > 0.69 in model 'owls.brm5'. We
recommend to run more iterations to get at least about 2200 posterior draws to
improve LOO-CV approximation accuracy.

Computed from 1500 by 599 log-likelihood matrix.

         Estimate    SE
elpd_loo  -1931.7  53.2
p_loo       185.2  13.2
looic      3863.3 106.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.69]   (good)     573   95.7%   83      
   (0.69, 1]   (bad)       25    4.2%   <NA>    
    (1, Inf)   (very bad)   1    0.2%   <NA>    
See help('pareto-k-diagnostic') for details.

(l.3 <- loo::loo(owls.brm6))

Warning: Found 29 observations with a pareto_k > 0.65 in model 'owls.brm6'. We
recommend to run more iterations to get at least about 2200 posterior draws to
improve LOO-CV approximation accuracy.

Computed from 750 by 599 log-likelihood matrix.

         Estimate    SE
elpd_loo  -1919.4  53.8
p_loo       183.6  12.7
looic      3838.9 107.6
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.65]   (good)     570   95.2%   32      
   (0.65, 1]   (bad)       27    4.5%   <NA>    
    (1, Inf)   (very bad)   2    0.3%   <NA>    
See help('pareto-k-diagnostic') for details.

(l.4 <- loo::loo(owls.brm7))

Warning: Found 9 observations with a pareto_k > 0.65 in model 'owls.brm7'. We
recommend to run more iterations to get at least about 2200 posterior draws to
improve LOO-CV approximation accuracy.

Computed from 750 by 599 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1665.1 28.3
p_loo        50.3  4.7
looic      3330.1 56.5
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.3]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.65]   (good)     590   98.5%   96      
   (0.65, 1]   (bad)        8    1.3%   <NA>    
    (1, Inf)   (very bad)   1    0.2%   <NA>    
See help('pareto-k-diagnostic') for details.

(l.5 <- loo::loo(owls.brm8))

Warning: Found 6 observations with a pareto_k > 0.7 in model 'owls.brm8'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.

Computed from 1500 by 599 log-likelihood matrix.

         Estimate   SE
elpd_loo  -1656.5 29.0
p_loo        50.1  4.0
looic      3312.9 57.9
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).

Pareto k diagnostic values:
                          Count Pct.    Min. ESS
(-Inf, 0.69]   (good)     593   99.0%   241     
   (0.69, 1]   (bad)        6    1.0%   <NA>    
    (1, Inf)   (very bad)   0    0.0%   <NA>    
See help('pareto-k-diagnostic') for details.

loo_compare(l.1, l.2, l.3, l.4, l.5)

          elpd_diff se_diff
owls.brm8    0.0       0.0 
owls.brm7   -8.6       5.1 
owls.brm6 -263.0      34.8 
owls.brm5 -275.2      34.9 
owls.brm4 -770.6      65.7 

10 Partial effects plots

brms

conditional_effects

owls.brm8 |>
  conditional_effects("FoodTreatment:SexParent") |>
  ## plot(points = TRUE, point.args = list(width=0.25))
  ## plot(points = TRUE, jitter_width = c(0.25,0))
  plot(points = TRUE, jitter_width = 0.25)

Warning: 'jitter_width' is deprecated. Please use 'point_args = list(width =
<width>)' instead.

These predictions appear to be based on the mean BroodSize of approximately 4.39. That is, the predictions are for a nest, not per chick. There does not appear to be a way to indicate the offset value.

ggpredict

owls.brm8 |>
  ggpredict(~ FoodTreatment * SexParent) |>
  plot(show_data = TRUE, jitter = c(0.25, 0))

Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.

Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail

These predictions appear to be based on the mean BroodSize of approximately 4.39. That is, the predictions are for a nest, not per chick. There does not appear to be a way to indicate the offset value.

ggemmeans

ggemmeans() can accommodate the offset correctly. There are two sensible choices:

• set the offset to the (log) of the mean BroodSize (similar to other partial effects), hence giving predictions appropriate for the average brood size conclusion.

off <- owls |> summarize(Mean = mean(BroodSize))
## as.numeric(off)
## owls.brm7 |>
##     ggemmeans(~FoodTreatment*SexParent, offset=off$Mean) |>
##     plot()
## owls.brm7 |>
##     ggemmeans(~FoodTreatment*SexParent, offset=log(off$Mean)) |>
##     plot()
## owls.brm8 |>
##     ggemmeans(~FoodTreatment*SexParent, offset=log(1)) |>
##     plot()
owls.brm8 |>
  ggemmeans(~ FoodTreatment * SexParent, offset = log(off$Mean)) |>
  plot()

Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, ziformula = ziformula, :
some components missing from 'family': downstream methods may fail

11 Model investigation

brms

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

owls.brm8 |> summary()

 Family: zero_inflated_negbinomial 
  Links: mu = log; shape = identity; zi = logit 
Formula: NCalls ~ FoodTreatment * SexParent + offset(log(BroodSize)) + (FoodTreatment * SexParent | Nest) 
         zi ~ FoodTreatment * SexParent
   Data: owls (Number of observations: 599) 
  Draws: 3 chains, each with iter = 10000; warmup = 5000; thin = 10;
         total post-warmup draws = 1500

Multilevel Hyperparameters:
~Nest (Number of levels: 27) 
                                                               Estimate
sd(Intercept)                                                      0.30
sd(FoodTreatmentSatiated)                                          1.11
sd(SexParentMale)                                                  0.17
sd(FoodTreatmentSatiated:SexParentMale)                            0.35
cor(Intercept,FoodTreatmentSatiated)                               0.09
cor(Intercept,SexParentMale)                                      -0.17
cor(FoodTreatmentSatiated,SexParentMale)                          -0.22
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                -0.11
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)    -0.04
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)            -0.06
                                                               Est.Error
sd(Intercept)                                                       0.10
sd(FoodTreatmentSatiated)                                           0.27
sd(SexParentMale)                                                   0.12
sd(FoodTreatmentSatiated:SexParentMale)                             0.22
cor(Intercept,FoodTreatmentSatiated)                                0.28
cor(Intercept,SexParentMale)                                        0.42
cor(FoodTreatmentSatiated,SexParentMale)                            0.40
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                  0.41
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)      0.43
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)              0.44
                                                               l-95% CI
sd(Intercept)                                                      0.13
sd(FoodTreatmentSatiated)                                          0.62
sd(SexParentMale)                                                  0.01
sd(FoodTreatmentSatiated:SexParentMale)                            0.02
cor(Intercept,FoodTreatmentSatiated)                              -0.45
cor(Intercept,SexParentMale)                                      -0.84
cor(FoodTreatmentSatiated,SexParentMale)                          -0.88
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                -0.83
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)    -0.81
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)            -0.82
                                                               u-95% CI Rhat
sd(Intercept)                                                      0.51 1.00
sd(FoodTreatmentSatiated)                                          1.72 1.00
sd(SexParentMale)                                                  0.43 1.00
sd(FoodTreatmentSatiated:SexParentMale)                            0.85 1.00
cor(Intercept,FoodTreatmentSatiated)                               0.63 1.00
cor(Intercept,SexParentMale)                                       0.72 1.00
cor(FoodTreatmentSatiated,SexParentMale)                           0.63 1.00
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                 0.70 1.00
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)     0.79 1.00
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)             0.78 1.00
                                                               Bulk_ESS
sd(Intercept)                                                      1462
sd(FoodTreatmentSatiated)                                          1560
sd(SexParentMale)                                                  1311
sd(FoodTreatmentSatiated:SexParentMale)                            1534
cor(Intercept,FoodTreatmentSatiated)                               1321
cor(Intercept,SexParentMale)                                       1508
cor(FoodTreatmentSatiated,SexParentMale)                           1323
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                 1486
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)     1570
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)             1446
                                                               Tail_ESS
sd(Intercept)                                                      1160
sd(FoodTreatmentSatiated)                                          1543
sd(SexParentMale)                                                  1522
sd(FoodTreatmentSatiated:SexParentMale)                            1342
cor(Intercept,FoodTreatmentSatiated)                               1260
cor(Intercept,SexParentMale)                                       1545
cor(FoodTreatmentSatiated,SexParentMale)                           1420
cor(Intercept,FoodTreatmentSatiated:SexParentMale)                 1276
cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale)     1451
cor(SexParentMale,FoodTreatmentSatiated:SexParentMale)             1454

Regression Coefficients:
                                       Estimate Est.Error l-95% CI u-95% CI
Intercept                                  0.83      0.10     0.64     1.03
zi_Intercept                              -1.60      0.25    -2.12    -1.15
FoodTreatmentSatiated                     -0.74      0.28    -1.31    -0.20
SexParentMale                             -0.09      0.11    -0.32     0.12
FoodTreatmentSatiated:SexParentMale        0.11      0.21    -0.32     0.51
zi_FoodTreatmentSatiated                   0.82      0.34     0.17     1.48
zi_SexParentMale                          -0.80      0.36    -1.52    -0.11
zi_FoodTreatmentSatiated:SexParentMale     0.61      0.43    -0.22     1.49
                                       Rhat Bulk_ESS Tail_ESS
Intercept                              1.00     1552     1523
zi_Intercept                           1.00     1238     1412
FoodTreatmentSatiated                  1.00     1409     1442
SexParentMale                          1.00     1532     1609
FoodTreatmentSatiated:SexParentMale    1.00     1462     1418
zi_FoodTreatmentSatiated               1.00     1258     1230
zi_SexParentMale                       1.00     1260     1379
zi_FoodTreatmentSatiated:SexParentMale 1.00     1388     1383

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     2.73      0.33     2.14     3.44 1.00     1308     1356

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• the average number of calls made per food deprived owl chicks when the female parent returns to the nest is 0.83 (on the link scale). When back-transformed, this is 2.29.
• on average, satiated chicks negotiate -1 (link scale) less when the female returns than deprived chicks. This equates to 0.48 fold fewer times and represents a 52% decline.
• on average, when the male parent returns, deprived chicks negotiate 0 (link scale) less when the female returns. This equates to 0.91 fold fewer times and represents a 9% decline (although this is not significant).
• there was not significantly detectable interaction suggesting that the effect of food treatment was consistent across both parents and vice versa.
• the estimated rate of non-detection of calls (false zeros) for deprived chicks when the female parent returns is -1.6 (logit scale). When back-transformed to the odds scale, this equates to the odds of a zero being false are 0.2:1. Expressed on a probability scale, it would suggest that the probability of a zero being false is 0.17.
• when the chicks are satiated, the odds not detecting a call (false zero) are increased 2.26 fold. That is, the odds increase by 30.65%.
• when the returning parent is male (for deprived chicks), the odds of not detecting a call (false zeros) are reduced 0.45 fold. That is, the odds decline by 68.96%.
• there is substantially more variation in sibling negotiations between food treatment effects within a nest than either between nests or between the sex of the parent effects within nests.
• variation in the sex of the parent effect is slightly negatively correlated to the variation between nests
• variation in the sex of the parent effect is negatively correlated to the variation in the interaction effect.

as_draws_df (posteriors)

# A draws_df: 500 iterations, 3 chains, and 138 variables
   b_Intercept b_zi_Intercept b_FoodTreatmentSatiated b_SexParentMale
1         0.68           -2.1                   -0.49           0.201
2         0.96           -1.9                   -1.04          -0.312
3         0.58           -1.7                   -0.69           0.189
4         0.83           -1.1                   -1.19          -0.073
5         1.04           -1.6                   -0.53          -0.190
6         0.69           -1.4                   -0.61          -0.097
7         0.72           -1.8                   -1.01           0.043
8         0.94           -1.9                   -0.93          -0.238
9         0.98           -1.7                   -1.07          -0.256
10        0.83           -1.7                   -0.61          -0.119
   b_FoodTreatmentSatiated:SexParentMale b_zi_FoodTreatmentSatiated
1                                  -0.82                       1.23
2                                   0.45                       1.15
3                                   0.14                       0.83
4                                   0.19                       0.45
5                                   0.22                       1.04
6                                   0.19                       0.47
7                                   0.22                       1.08
8                                   0.21                       1.15
9                                   0.39                       0.81
10                                 -0.03                       0.89
   b_zi_SexParentMale b_zi_FoodTreatmentSatiated:SexParentMale
1               -0.59                                    0.511
2               -0.35                                   -0.096
3               -0.78                                    0.748
4               -0.84                                    0.623
5               -0.62                                    0.481
6               -0.85                                    0.623
7               -0.46                                    0.348
8               -0.30                                    0.107
9               -0.87                                    0.955
10              -0.11                                   -0.321
# ... with 1490 more draws, and 130 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

# A tibble: 138 × 10
   variable     median lower upper      Pl     Pg  rhat length ess_bulk ess_tail
   <chr>         <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>  <dbl>    <dbl>    <dbl>
 1 b_Intercept   2.28  1.89  2.80  0       1      1.00    1500    1551.    1523.
 2 b_zi_Interc…  0.206 0.117 0.309 1       0      1.00    1500    1236.    1412.
 3 b_FoodTreat…  0.477 0.231 0.748 0.996   0.004  1.00    1500    1409.    1442.
 4 b_SexParent…  0.915 0.725 1.12  0.805   0.195  1.00    1500    1532.    1609.
 5 b_FoodTreat…  1.11  0.697 1.61  0.281   0.719  1.00    1500    1463.    1418.
 6 b_zi_FoodTr…  2.25  1.04  4.10  0.00933 0.991  1.000   1500    1258.    1230.
 7 b_zi_SexPar…  0.456 0.186 0.821 0.987   0.0133 1.00    1500    1260.    1379.
 8 b_zi_FoodTr…  1.82  0.566 3.83  0.07    0.93   1.00    1500    1388.    1383.
 9 sd_Nest__In…  1.35  1.13  1.64  0       1      1.00    1500    1461.    1160.
10 sd_Nest__Fo…  2.95  1.76  5.20  0       1      1.00    1500    1560.    1543.
# ℹ 128 more rows

[1] 0.1701245

[1] 0.3166131

[1] 0.08428682

~17% of zeros are false zeros = 1/detectability

tidyMCMC

owls.brm8$fit |>
  tidyMCMC(
    estimate.method = "median",
    conf.int = TRUE, conf.method = "HPDinterval",
    rhat = TRUE, ess = TRUE
  )

# A tibble: 137 × 7
   term                        estimate std.error conf.low conf.high  rhat   ess
   <chr>                          <dbl>     <dbl>    <dbl>     <dbl> <dbl> <int>
 1 b_Intercept                   0.827     0.101     0.637     1.03  1.00   1558
 2 b_zi_Intercept               -1.60      0.254    -2.13     -1.16  1.000  1213
 3 b_FoodTreatmentSatiated      -0.740     0.282    -1.31     -0.221 1.00   1424
 4 b_SexParentMale              -0.0926    0.110    -0.292     0.138 1.00   1535
 5 b_FoodTreatmentSatiated:Se…   0.107     0.205    -0.294     0.530 1.00   1463
 6 b_zi_FoodTreatmentSatiated    0.816     0.335     0.165     1.48  0.999  1209
 7 b_zi_SexParentMale           -0.798     0.357    -1.57     -0.171 1.00   1274
 8 b_zi_FoodTreatmentSatiated…   0.611     0.426    -0.194     1.51  0.999  1374
 9 sd_Nest__Intercept            0.303     0.0958    0.126     0.498 1.000  1457
10 sd_Nest__FoodTreatmentSati…   1.11      0.275     0.589     1.67  1.00   1543
# ℹ 127 more rows

Conclusions:

see above

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

owls.brm8 |> get_variables()

  [1] "b_Intercept"                                                         
  [2] "b_zi_Intercept"                                                      
  [3] "b_FoodTreatmentSatiated"                                             
  [4] "b_SexParentMale"                                                     
  [5] "b_FoodTreatmentSatiated:SexParentMale"                               
  [6] "b_zi_FoodTreatmentSatiated"                                          
  [7] "b_zi_SexParentMale"                                                  
  [8] "b_zi_FoodTreatmentSatiated:SexParentMale"                            
  [9] "sd_Nest__Intercept"                                                  
 [10] "sd_Nest__FoodTreatmentSatiated"                                      
 [11] "sd_Nest__SexParentMale"                                              
 [12] "sd_Nest__FoodTreatmentSatiated:SexParentMale"                        
 [13] "cor_Nest__Intercept__FoodTreatmentSatiated"                          
 [14] "cor_Nest__Intercept__SexParentMale"                                  
 [15] "cor_Nest__FoodTreatmentSatiated__SexParentMale"                      
 [16] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"            
 [17] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
 [18] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"        
 [19] "shape"                                                               
 [20] "Intercept"                                                           
 [21] "Intercept_zi"                                                        
 [22] "r_Nest[AutavauxTV,Intercept]"                                        
 [23] "r_Nest[Bochet,Intercept]"                                            
 [24] "r_Nest[Champmartin,Intercept]"                                       
 [25] "r_Nest[ChEsard,Intercept]"                                           
 [26] "r_Nest[Chevroux,Intercept]"                                          
 [27] "r_Nest[CorcellesFavres,Intercept]"                                   
 [28] "r_Nest[Etrabloz,Intercept]"                                          
 [29] "r_Nest[Forel,Intercept]"                                             
 [30] "r_Nest[Franex,Intercept]"                                            
 [31] "r_Nest[GDLV,Intercept]"                                              
 [32] "r_Nest[Gletterens,Intercept]"                                        
 [33] "r_Nest[Henniez,Intercept]"                                           
 [34] "r_Nest[Jeuss,Intercept]"                                             
 [35] "r_Nest[LesPlanches,Intercept]"                                       
 [36] "r_Nest[Lucens,Intercept]"                                            
 [37] "r_Nest[Lully,Intercept]"                                             
 [38] "r_Nest[Marnand,Intercept]"                                           
 [39] "r_Nest[Montet,Intercept]"                                            
 [40] "r_Nest[Murist,Intercept]"                                            
 [41] "r_Nest[Oleyes,Intercept]"                                            
 [42] "r_Nest[Payerne,Intercept]"                                           
 [43] "r_Nest[Rueyes,Intercept]"                                            
 [44] "r_Nest[Seiry,Intercept]"                                             
 [45] "r_Nest[Sevaz,Intercept]"                                             
 [46] "r_Nest[StAubin,Intercept]"                                           
 [47] "r_Nest[Trey,Intercept]"                                              
 [48] "r_Nest[Yvonnand,Intercept]"                                          
 [49] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"                            
 [50] "r_Nest[Bochet,FoodTreatmentSatiated]"                                
 [51] "r_Nest[Champmartin,FoodTreatmentSatiated]"                           
 [52] "r_Nest[ChEsard,FoodTreatmentSatiated]"                               
 [53] "r_Nest[Chevroux,FoodTreatmentSatiated]"                              
 [54] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"                       
 [55] "r_Nest[Etrabloz,FoodTreatmentSatiated]"                              
 [56] "r_Nest[Forel,FoodTreatmentSatiated]"                                 
 [57] "r_Nest[Franex,FoodTreatmentSatiated]"                                
 [58] "r_Nest[GDLV,FoodTreatmentSatiated]"                                  
 [59] "r_Nest[Gletterens,FoodTreatmentSatiated]"                            
 [60] "r_Nest[Henniez,FoodTreatmentSatiated]"                               
 [61] "r_Nest[Jeuss,FoodTreatmentSatiated]"                                 
 [62] "r_Nest[LesPlanches,FoodTreatmentSatiated]"                           
 [63] "r_Nest[Lucens,FoodTreatmentSatiated]"                                
 [64] "r_Nest[Lully,FoodTreatmentSatiated]"                                 
 [65] "r_Nest[Marnand,FoodTreatmentSatiated]"                               
 [66] "r_Nest[Montet,FoodTreatmentSatiated]"                                
 [67] "r_Nest[Murist,FoodTreatmentSatiated]"                                
 [68] "r_Nest[Oleyes,FoodTreatmentSatiated]"                                
 [69] "r_Nest[Payerne,FoodTreatmentSatiated]"                               
 [70] "r_Nest[Rueyes,FoodTreatmentSatiated]"                                
 [71] "r_Nest[Seiry,FoodTreatmentSatiated]"                                 
 [72] "r_Nest[Sevaz,FoodTreatmentSatiated]"                                 
 [73] "r_Nest[StAubin,FoodTreatmentSatiated]"                               
 [74] "r_Nest[Trey,FoodTreatmentSatiated]"                                  
 [75] "r_Nest[Yvonnand,FoodTreatmentSatiated]"                              
 [76] "r_Nest[AutavauxTV,SexParentMale]"                                    
 [77] "r_Nest[Bochet,SexParentMale]"                                        
 [78] "r_Nest[Champmartin,SexParentMale]"                                   
 [79] "r_Nest[ChEsard,SexParentMale]"                                       
 [80] "r_Nest[Chevroux,SexParentMale]"                                      
 [81] "r_Nest[CorcellesFavres,SexParentMale]"                               
 [82] "r_Nest[Etrabloz,SexParentMale]"                                      
 [83] "r_Nest[Forel,SexParentMale]"                                         
 [84] "r_Nest[Franex,SexParentMale]"                                        
 [85] "r_Nest[GDLV,SexParentMale]"                                          
 [86] "r_Nest[Gletterens,SexParentMale]"                                    
 [87] "r_Nest[Henniez,SexParentMale]"                                       
 [88] "r_Nest[Jeuss,SexParentMale]"                                         
 [89] "r_Nest[LesPlanches,SexParentMale]"                                   
 [90] "r_Nest[Lucens,SexParentMale]"                                        
 [91] "r_Nest[Lully,SexParentMale]"                                         
 [92] "r_Nest[Marnand,SexParentMale]"                                       
 [93] "r_Nest[Montet,SexParentMale]"                                        
 [94] "r_Nest[Murist,SexParentMale]"                                        
 [95] "r_Nest[Oleyes,SexParentMale]"                                        
 [96] "r_Nest[Payerne,SexParentMale]"                                       
 [97] "r_Nest[Rueyes,SexParentMale]"                                        
 [98] "r_Nest[Seiry,SexParentMale]"                                         
 [99] "r_Nest[Sevaz,SexParentMale]"                                         
[100] "r_Nest[StAubin,SexParentMale]"                                       
[101] "r_Nest[Trey,SexParentMale]"                                          
[102] "r_Nest[Yvonnand,SexParentMale]"                                      
[103] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"              
[104] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"                  
[105] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"             
[106] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"                 
[107] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"                
[108] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"         
[109] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"                
[110] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"                   
[111] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"                  
[112] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"                    
[113] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"              
[114] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"                 
[115] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"                   
[116] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"             
[117] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"                  
[118] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"                   
[119] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"                 
[120] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"                  
[121] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"                  
[122] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"                  
[123] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"                 
[124] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"                  
[125] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"                   
[126] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"                   
[127] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"                 
[128] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"                    
[129] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"                
[130] "prior_Intercept"                                                     
[131] "prior_b"                                                             
[132] "prior_shape"                                                         
[133] "prior_b_zi"                                                          
[134] "prior_Intercept_zi"                                                  
[135] "prior_sd_Nest"                                                       
[136] "prior_cor_Nest"                                                      
[137] "lprior"                                                              
[138] "lp__"                                                                
[139] "accept_stat__"                                                       
[140] "stepsize__"                                                          
[141] "treedepth__"                                                         
[142] "n_leapfrog__"                                                        
[143] "divergent__"                                                         
[144] "energy__"                                                            

owls.draw <- owls.brm8 |>
  gather_draws(`b.Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE)
owls.draw

# A tibble: 6,000 × 5
# Groups:   .variable [4]
   .chain .iteration .draw .variable   .value
    <int>      <int> <int> <chr>        <dbl>
 1      1          1     1 b_Intercept  0.682
 2      1          2     2 b_Intercept  0.956
 3      1          3     3 b_Intercept  0.577
 4      1          4     4 b_Intercept  0.827
 5      1          5     5 b_Intercept  1.04 
 6      1          6     6 b_Intercept  0.689
 7      1          7     7 b_Intercept  0.724
 8      1          8     8 b_Intercept  0.941
 9      1          9     9 b_Intercept  0.981
10      1         10    10 b_Intercept  0.834
# ℹ 5,990 more rows

We can then summarise this

owls.draw |> median_hdci()

# A tibble: 4 × 7
  .variable                         .value .lower .upper .width .point .interval
  <chr>                              <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_FoodTreatmentSatiated          -0.741  -1.31  -0.221   0.95 median hdci     
2 b_FoodTreatmentSatiated:SexPare…  0.104  -0.294  0.530   0.95 median hdci     
3 b_Intercept                       0.824   0.637  1.03    0.95 median hdci     
4 b_SexParentMale                  -0.0891 -0.304  0.127   0.95 median hdci     

## On a fractional scale
owls.draw |>
  mutate(.value = exp(.value)) |>
  median_hdci()

# A tibble: 4 × 7
  .variable                         .value .lower .upper .width .point .interval
  <chr>                              <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 b_FoodTreatmentSatiated            0.477  0.245  0.770   0.95 median hdci     
2 b_FoodTreatmentSatiated:SexParen…  1.11   0.697  1.61    0.95 median hdci     
3 b_Intercept                        2.28   1.89   2.80    0.95 median hdci     
4 b_SexParentMale                    0.915  0.738  1.14    0.95 median hdci     

Lets plot the parameter posteriors (on the link scale - since we are including the intercept).

owls.brm8 |>
  gather_draws(`b_Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE) |>
  ## mutate(.value = exp(.value)) |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_slab(aes(
    x = .value, y = .variable,
    fill = stat(ggdist::cut_cdf_qi(cdf,
      .width = c(0.5, 0.8, 0.95),
      labels = scales::percent_format()
    ))
  ), color = "black") +
  scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)

Warning: `stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95), labels =
scales::percent_format()))` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95),
  labels = scales::percent_format()))` instead.

owls.brm8 |>
  gather_draws(`.Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE) |>
  ggplot() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  stat_halfeye(aes(x = .value, y = .variable)) +
  theme_classic()

bayesplot

owls.brm8$fit |> plot(type = "intervals")

'pars' not specified. Showing first 10 parameters by default.

ci_level: 0.8 (80% intervals)

outer_level: 0.95 (95% intervals)

half-eye (ggdist)

owls.brm8 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  facet_wrap(~.variable, scales = "free") +
  theme(axis.text.y = element_blank())

owls.brm8 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(.variable != "b_Intercept") |>
  ggplot() +
  stat_halfeye(aes(x = .value, y = .variable)) +
  geom_vline(xintercept = 0, linetype = "dashed")

density ridges (ggridges)

owls.brm8 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(str_detect(.variable, "b_.*Intercept", negate = TRUE)) |>
  ggplot() +
  geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed")

Picking joint bandwidth of 0.0588

## Or in colour
owls.brm8 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(str_detect(.variable, "b_.*Intercept", negate = TRUE)) |>
  ggplot() +
  geom_density_ridges_gradient(
    aes(
      x = (.value),
      y = .variable,
      fill = stat(x)
    ),
    alpha = 0.4, colour = "white",
    quantile_lines = TRUE,
    quantiles = c(0.025, 0.975)
  ) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous() +
  scale_fill_viridis_c(option = "C")

Picking joint bandwidth of 0.0588

## Fractional scale
owls.brm8 |>
  gather_draws(`^b_.*`, regex = TRUE) |>
  filter(str_detect(.variable, "b_.*Intercept", negate = TRUE)) |>
  ggplot() +
  geom_density_ridges_gradient(
    aes(
      x = exp(.value),
      y = .variable,
      fill = stat(x)
    ),
    alpha = 0.4, colour = "white",
    quantile_lines = TRUE,
    quantiles = c(0.025, 0.975)
  ) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  scale_x_continuous(trans = scales::log2_trans()) +
  scale_fill_viridis_c(option = "C")

Picking joint bandwidth of 0.0848

tidy_draws

This is purely a graphical depiction on the posteriors.

owls.brm8 |> tidy_draws()

# A tibble: 1,500 × 147
   .chain .iteration .draw b_Intercept b_zi_Intercept b_FoodTreatmentSatiated
    <int>      <int> <int>       <dbl>          <dbl>                   <dbl>
 1      1          1     1       0.682          -2.07                  -0.487
 2      1          2     2       0.956          -1.89                  -1.04 
 3      1          3     3       0.577          -1.67                  -0.688
 4      1          4     4       0.827          -1.11                  -1.19 
 5      1          5     5       1.04           -1.57                  -0.526
 6      1          6     6       0.689          -1.38                  -0.607
 7      1          7     7       0.724          -1.80                  -1.01 
 8      1          8     8       0.941          -1.92                  -0.932
 9      1          9     9       0.981          -1.73                  -1.07 
10      1         10    10       0.834          -1.74                  -0.613
# ℹ 1,490 more rows
# ℹ 141 more variables: b_SexParentMale <dbl>,
#   `b_FoodTreatmentSatiated:SexParentMale` <dbl>,
#   b_zi_FoodTreatmentSatiated <dbl>, b_zi_SexParentMale <dbl>,
#   `b_zi_FoodTreatmentSatiated:SexParentMale` <dbl>, sd_Nest__Intercept <dbl>,
#   sd_Nest__FoodTreatmentSatiated <dbl>, sd_Nest__SexParentMale <dbl>,
#   `sd_Nest__FoodTreatmentSatiated:SexParentMale` <dbl>, …

spread_draws

owls.brm8 |> spread_draws(`.*Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE)

# A tibble: 1,500 × 43
   .chain .iteration .draw b_Intercept b_zi_Intercept b_FoodTreatmentSatiated
    <int>      <int> <int>       <dbl>          <dbl>                   <dbl>
 1      1          1     1       0.682          -2.07                  -0.487
 2      1          2     2       0.956          -1.89                  -1.04 
 3      1          3     3       0.577          -1.67                  -0.688
 4      1          4     4       0.827          -1.11                  -1.19 
 5      1          5     5       1.04           -1.57                  -0.526
 6      1          6     6       0.689          -1.38                  -0.607
 7      1          7     7       0.724          -1.80                  -1.01 
 8      1          8     8       0.941          -1.92                  -0.932
 9      1          9     9       0.981          -1.73                  -1.07 
10      1         10    10       0.834          -1.74                  -0.613
# ℹ 1,490 more rows
# ℹ 37 more variables: b_SexParentMale <dbl>,
#   `b_FoodTreatmentSatiated:SexParentMale` <dbl>, sd_Nest__Intercept <dbl>,
#   cor_Nest__Intercept__FoodTreatmentSatiated <dbl>,
#   cor_Nest__Intercept__SexParentMale <dbl>,
#   `cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale` <dbl>,
#   Intercept <dbl>, Intercept_zi <dbl>, …

posterior_samples

owls.brm8 |>
  posterior_samples() |>
  as_tibble()

Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.

# A tibble: 1,500 × 138
   b_Intercept b_zi_Intercept b_FoodTreatmentSatiated b_SexParentMale
         <dbl>          <dbl>                   <dbl>           <dbl>
 1       0.682          -2.07                  -0.487          0.201 
 2       0.956          -1.89                  -1.04          -0.312 
 3       0.577          -1.67                  -0.688          0.189 
 4       0.827          -1.11                  -1.19          -0.0725
 5       1.04           -1.57                  -0.526         -0.190 
 6       0.689          -1.38                  -0.607         -0.0969
 7       0.724          -1.80                  -1.01           0.0427
 8       0.941          -1.92                  -0.932         -0.238 
 9       0.981          -1.73                  -1.07          -0.256 
10       0.834          -1.74                  -0.613         -0.119 
# ℹ 1,490 more rows
# ℹ 134 more variables: `b_FoodTreatmentSatiated:SexParentMale` <dbl>,
#   b_zi_FoodTreatmentSatiated <dbl>, b_zi_SexParentMale <dbl>,
#   `b_zi_FoodTreatmentSatiated:SexParentMale` <dbl>, sd_Nest__Intercept <dbl>,
#   sd_Nest__FoodTreatmentSatiated <dbl>, sd_Nest__SexParentMale <dbl>,
#   `sd_Nest__FoodTreatmentSatiated:SexParentMale` <dbl>,
#   cor_Nest__Intercept__FoodTreatmentSatiated <dbl>, …

\(R^2\)

owls.brm8 |>
  bayes_R2(re.form = NA, summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.1732165 0.1088131 0.2455375   0.95 median      hdci

owls.brm8 |>
  bayes_R2(re.form = ~ (1 | Nest), summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.2024367 0.1311129 0.2833968   0.95 median      hdci

owls.brm8 |>
  bayes_R2(re.form = ~ (FoodTreatment * SexParent | Nest), summary = FALSE) |>
  median_hdci()

          y      ymin      ymax .width .point .interval
1 0.2419584 0.1894629 0.3002374   0.95 median      hdci

Modelsummary

owls.brm8 |> modelsummary(
  statistic = c("conf.low", "conf.high"),
  shape = term ~ statistic,
  exponentiate = TRUE
)

owls.brm8 |> modelplot(exponentiate = TRUE)

12 Further investigations

brms

## ## The following should work, but there is a bug and therefore it does not
## ## (although it has been reported - so may get fixed at some point).
## ## The offset seems to get handled incorrectly
## newdata <- owls.brm8 |>
##     emmeans(~FoodTreatment|SexParent, offset=0, type='response') |>
##     as.data.frame()
newdata <- owls.brm8 |>
  emmeans(~ FoodTreatment | SexParent, type = "response") |>
  as.data.frame()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

head(newdata)

 FoodTreatment SexParent     prob lower.HPD upper.HPD
 Deprived      Female    2.278798 1.8900365  2.801815
 Satiated      Female    1.093940 0.5790027  1.752645
 Deprived      Male      2.085553 1.7539866  2.466215
 Satiated      Male      1.118972 0.6283216  1.754491

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

## As an alternative, we can do the following...
newdata <- owls.brm8 |>
  emmeans(~ FoodTreatment | SexParent,
    at = list(BroodSize = 1), type = "response"
  ) |>
  as.data.frame()

Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.

head(newdata)

 FoodTreatment SexParent     prob lower.HPD upper.HPD
 Deprived      Female    2.278798 1.8900365  2.801815
 Satiated      Female    1.093940 0.5790027  1.752645
 Deprived      Male      2.085553 1.7539866  2.466215
 Satiated      Male      1.118972 0.6283216  1.754491

Point estimate displayed: median 
Results are back-transformed from the log scale 
HPD interval probability: 0.95 

ggplot(newdata) +
  geom_pointrange(
    aes(
      y = prob, x = FoodTreatment, color = SexParent,
      ymin = lower.HPD, ymax = upper.HPD
    ),
    position = position_dodge(width = 0.2)
  ) +
  theme_classic()

12.1 rstanarm

12.2 rstanarm

12.3 brms

:::

13 References